1.1 The mixtures model

PCP models an observed exposure matrix \(D\) as the sum of three underlying ground-truth matrices: \[ \underset{\text{mixture}}{D_{n \times p}} = \underset{\text{low-rank}}{L_0} + \underset{\text{sparse}}{S_0} + \underset{\text{noise}}{Z_0} \] a low-rank matrix \(L_0\) encoding consistent patterns of exposure, a sparse matrix \(S_0\) isolating unique or outlying exposure events (that cannot be explained by the consistent exposure patterns), and dense noise \(Z_0\). All of these matrices are of dimension \(n \times p\), where \(n\) is the number of observations (e.g., study participants or measurement dates) and \(p\) is the number of exposures (chemical and/or non-chemical stressors). Beyond this mixtures model, the main (EH-relevant) assumption made by PCP is that \(Z_0 \sim \mathcal{N}(\mu, \sigma^2)\) consists of i.i.d. Gaussian noise corrupting each entry of the overall exposure matrix \(D\). A visual example of PCP’s mixtures model is provided below:

A visual example of the PCP mixtures model.

Under these minimal assumptions, PCP aims to decompose the observed mixture \(D\) into the above mixtures model via an optimization program, thereby recovering the ground-truth \(L_0\) and \(S_0\) with accurate estimates \(\hat{L}\) and \(\hat{S}\).

1.2 The low-rank matrix

The estimated low-rank matrix \(\hat{L}\) provides information on the consistent exposure patterns, satisfying \[ r = \text{rank}(\hat{L}) \ll \min(n, p). \] The rank \(r\) corresponds to the (relatively few) number of underlying patterns governing the mixture, such as specific sources or behaviors leading to exposure. Notice that \(\hat{L} \in \mathbb{R}^{n \times p}\), meaning it is still defined in terms of the original variables. Put differently, \(\hat{L}\) can be taken as a robust approximation to the true environmental mixture matrix, unperturbed by outliers (captured in \(\hat{S}\)) or noise (handled by PCP’s noise term). In this way, the latent exposure patterns are encoded in \(\hat{L}\) rather than directly estimated. To explicitly obtain the exposure patterns from \(\hat{L}\), PCP may then be paired with various matrix factorization methods (e.g., PCA, factor analysis, or NMF) that yield chemical loadings and individual scores for use in downstream health models.

1.3 The sparse matrix

The estimated sparse matrix \(\hat{S}\) captures unusually high/low outlying exposure events, unexplained by the identified patterns in \(\hat{L}\). Most entries in \(\hat{S}\) are 0, with non-zero entries identifying such extreme exposure activity. The number, location (i.e., support), and value of non-zero entries in \(\hat{S}\) need not be defined; PCP isolates these during optimization.

2.1 Root PCP

The convex model available in pcpr is Root PCP, denoted here as \(\sqrt{\text{PCP}}\). For a comprehensive technical understanding, we refer readers to Zhang et al. (2021) introducing \(\sqrt{\text{PCP}}\). Here we briefly introduce the root_pcp function in pcpr.

2.2 RRMC

3.1 Exploring the raw queens data

First, let’s load the pcpr package (along with a few other packages we will need for data processing) and take a look at the queens dataset:

library(dplyr)
library(ggplot2)
library(magrittr)
library(pcpr)
library(stringr)
library(tidyr)

queens 

## # A tibble: 2,443 × 27
##    Date            Al   NH4      As     Ba       Br     Cd      Ca      Cl
##    <date>       <dbl> <dbl>   <dbl>  <dbl>    <dbl>  <dbl>   <dbl>   <dbl>
##  1 2001-04-04 NA      1.62  NA      NA     NA       NA     NA      NA     
##  2 2001-04-07  0      2.66   0       0.012  0.00488  0      0.0401  0.0079
##  3 2001-04-13  0.0094 1.41   0.0016  0.024  0.00211  0.004  0.036   0     
##  4 2001-04-19  0.0104 1.22   0.001   0.006  0.00422  0      0.0543  0.003 
##  5 2001-04-25  0.0172 0.723  0.0024  0.015  0.00117  0      0.0398  0     
##  6 2001-05-01  0.0384 3.48   0.0017  0.041  0.00873  0.001  0.136   0     
##  7 2001-05-04  0.0964 6.22   0.0025  0.039  0.0111   0      0.137   0     
##  8 2001-05-07  0.004  0.233  0.001   0.016  0.00263  0      0.055   0.0054
##  9 2001-05-10  0.0547 2.04   0.001   0.055  0.00521  0      0.121   0.001 
## 10 2001-05-13  0.0215 0.229  0       0.021  0.00122  0      0.0249  0     
## # ℹ 2,433 more rows
## # ℹ 18 more variables: Cr <dbl>, Cu <dbl>, EC <dbl>, Fe <dbl>, Pb <dbl>,
## #   Mg <dbl>, Mn <dbl>, Ni <dbl>, OC <dbl>, K <dbl>, Se <dbl>, Si <dbl>,
## #   Na <dbl>, S <dbl>, Ti <dbl>, NO3 <dbl>, V <dbl>, Zn <dbl>

Let’s spend some time visualizing some of the trends of the raw dataset, to get a sense of what we’re dealing with. We can start by plotting each of the measured chemical species as timeseries. Plotted in black are the raw observed PM\(_{2.5}\) measurements (in µg/m\(^3\)) over time (04/04/2001 - 12/30/2021), and plotted in red are some (rough) lines of best fit:

queens %>%
  pivot_longer(
    colnames(queens)[-1], names_to = "chem", values_to = "concentration"
  ) %>%
  filter(!is.na(concentration)) %>% 
  ggplot(aes(x = Date, y = concentration)) +
  geom_line() +
  geom_smooth(color = "red", formula="y ~ x", method = "loess", span = 0.05) +
  facet_wrap(~chem, scales = "free_y") +
  labs(x = "Date", y = "Concentration (µg/m^3)") +
  theme_bw()

Some initial points to make note of, looking at the data:

• The y-axis scale for each chemical’s concentration differs widely (e.g., NH4 ranges from 0 to 9 µg/m\(^3\), while Ti ranges from 0 to 0.04 µg/m\(^3\)). We will have to take this into account during preprocessing (which we do in the next section Preliminary PCA).
• Some chemical species record negative measurements (e.g., Al, Ba, Cd). As explained in the EPA AQS data documentation, this is due to idiosyncrasies in the instruments used to collect the measurements.
• We can see strong seasonal trends in some species (e.g., NO3), as well as a reduction in concentration over time in others (e.g., Ni).
• Many species have prominent outliers, or extreme exposure events (e.g., Cr, Cu). Further, some species have extreme exposure events following seasonal trends, e.g., K.
• Elemental carbon (EC) and organic carbon (OC) are both missing measurements for the first 8 years of the dataset, likely because the air monitors for those two species were not operational from 2001 - 2009. Handling such systematic missingness is out of scope for this tutorial, so we make the simple decision to omit all data from 2001 - 2009, rather than dropping EC and OC, since including both will be helpful in our source apportionment study.

Below we remove the early measurement dates from 2001 to 2009 that are missing observations for EC and OC, yielding the queens_small dataset that we will use moving forward:

start_date_df <- queens %>% select(Date, EC, OC) %>% na.omit() %>% slice_head(n = 1)
start_date <- start_date_df$Date[1]
cat("Start date:", as.character(start_date))

## Start date: 2009-04-01

queens_small <- queens %>% filter(Date >= as.Date(start_date))
queens_small

## # A tibble: 1,535 × 27
##    Date          Al   NH4      As    Ba     Br    Cd     Ca     Cl    Cr     Cu
##    <date>     <dbl> <dbl>   <dbl> <dbl>  <dbl> <dbl>  <dbl>  <dbl> <dbl>  <dbl>
##  1 2009-04-01 0.044 0.883 0       0     0.0034 0     0.0507 0.677  0     0.0026
##  2 2009-04-04 0.008 0.899 0       0     0.0027 0.009 0.012  0      0.001 0.001 
##  3 2009-04-07 0.03  0.623 0       0     0      0     0.0142 0.006  0     0.0039
##  4 2009-04-10 0.016 3.43  0       0     0.0063 0     0.0674 0.0561 0     0.0062
##  5 2009-04-13 0.025 0.774 0.001   0     0.003  0.012 0.0372 0.001  0     0.0025
##  6 2009-04-16 0.002 0.822 0       0.003 0.0051 0.004 0.0458 0      0     0.001 
##  7 2009-04-19 0     0.748 0.00105 0     0.0015 0     0.0299 0.0547 0.002 0.001 
##  8 2009-04-22 0     2.1   0       0     0.0033 0     0.0488 0.0101 0.002 0.0051
##  9 2009-04-25 0.044 2.47  0       0     0.0037 0     0.0626 0.024  0.001 0.0064
## 10 2009-04-28 0.045 1.33  0       0     0.0043 0     0.0494 0.004  0     0.0023
## # ℹ 1,525 more rows
## # ℹ 16 more variables: EC <dbl>, Fe <dbl>, Pb <dbl>, Mg <dbl>, Mn <dbl>,
## #   Ni <dbl>, OC <dbl>, K <dbl>, Se <dbl>, Si <dbl>, Na <dbl>, S <dbl>,
## #   Ti <dbl>, NO3 <dbl>, V <dbl>, Zn <dbl>

Next, let’s take a look at the correlation structure of our queens_small data:

queens_small %>% 
  select(-Date) %>%
  as_tibble() %>% 
  GGally::ggcorr(., method = "pairwise.complete.obs", limits = F, label = F, size = 5)

The high dimensionality of our queens_small air pollution mixture matrix gives rise to this relatively complex correlation matrix. No strong patterns jump out right away here, although Na’s correlation with Cl and Mg calls to mind sea salt from the Atlantic Ocean or perhaps road salt kicked up by traffic. We’d like to employ PCP in order to reduce the complexity of our data for more robust downstream analysis. Keep this correlation matrix in mind, since after applying PCP, we’ll examine the correlations matrix of \(\hat{L}\) and compare.

3.2 Preliminary PCA

Before applying PCP to our mixture, it’s good to first take a look at what Principal Component Analysis (PCA) is able to extract from the queens_small data. We do this for a couple of reasons:

1. First, we’d like to establish a baseline for comparison. To better understand the effect of applying PCP to the queens_small data, we’d like to see what low dimensional structures exist in the raw data from the outset.
2. Second, understanding the behavior of the singular values governing the raw data will help inform choices we must make pertaining to PCP’s optimization scheme. A preliminary PCA will enable us to scrutinize our data’s singular values. \(\sqrt{\text{PCP}}\) is well-suited for data exhibiting rapidly decaying singular values (e.g., imaging data), while RRMC is best for more messy data with slowly decaying singular values (environmental mixtures data typically falls into this bucket).

Before we are able to call upon PCA (or PCP, for that matter) we need to first preprocess our data for better numerical stability. In practical terms, both PCA and PCP’s statistical routines are sensitive to data with variable scales. Before arbitrarily normalizing or standardizing our data, let’s examine the distribution of our queens_small dataset:

queens_small %>%
  pivot_longer(
    colnames(queens_small)[-1], names_to = "chem", values_to = "concentration"
  ) %>%
  filter(!is.na(concentration)) %>% 
  ggplot(aes(x = concentration)) +
  geom_histogram(bins = 50) +
  theme_bw() +
  facet_wrap(~chem, scales = "free")

Most of the chemical species look as though they are normally distributed. We have many choices for how we’d like to preprocess our data, e.g., standardize the data, min-max normalization, etc. Because most of our data appears (roughly) normally distributed, and because in PCA analyses, standardization is typical, let’s go ahead and standardize (scale and center) our data to have \(\mu = 0\), \(\sigma = 1\):

queens_standardized <- queens_small %>%
  select(-Date) %>%
  scale() %>%
  as_tibble() 

queens_standardized$Date <- queens_small$Date

queens_standardized %>% 
  pivot_longer(
    colnames(queens_small)[-1], names_to = "chem", values_to = "concentration"
  ) %>%
  filter(!is.na(concentration)) %>%
  ggplot(aes(x = concentration)) +
  geom_histogram(bins = 50) +
  theme_bw() +
  facet_wrap(~chem, scales = "free")

Great, now we’re ready to run the preliminary PCA analysis! We’ll use the stats::prcomp function for this, which by default omits missing values in our dataset. Instead, let’s handle NAs by imputing with the column-wise means. Below, we define two helper functions, mean_impute and pca:

mean_impute <- function(x) {
  x[is.na(x)] <- mean(x, na.rm = T)
  x
}

pca <- function(mat, pcs = c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6"), colgroups = NULL, impute = T) {
  
  if ("Date" %in% colnames(mat)) mat <- mat %>% select(-Date)
  if (impute) mat <- mat %>% apply(., 2, mean_impute)
  
  #### PCA ####
  pca.ln <- prcomp(mat) 
  
  #### VARIANCE TABLE ####
  singvalues_ln <- matrix(pca.ln$sdev^2)
  perc_variance <- round(100*matrix(pca.ln$sdev^2/sum(pca.ln$sdev^2)),1)
  
  pca_summary <- data.frame("Principle component" = 1:min(dim(mat)), "Singular values" = singvalues_ln, "Percent variance" = perc_variance, "Total cumulative variance" = purrr::accumulate(perc_variance, sum))
  
  var_tbl <- kableExtra::kbl(pca_summary, col.names = c("Principle component", "Singular value", "% variance", "Total cumulative variance"), align = "c") %>% 
    kableExtra::kable_classic(full_width = F, html_font = "Cambria", position = "center") %>% 
    kableExtra::kable_styling(bootstrap_options = c("hover", "condensed"), fixed_thead = T)
  
  #### LOADINGS ####
  pca.ln.ld <- as.data.frame.matrix(pca.ln$rotation)
  pca.ln.ld$chem <- row.names(pca.ln.ld)
  
  if (!is.null(colgroups)) {
    colgroups <- colgroups %>% dplyr::rename(chem = !!names(colgroups)[1])
  } else {
    colgroups <- data.frame(chem = colnames(mat), group = "1")
  }
  grouping <- names(colgroups)[2]
  
  plot_loadings_pca <- pca.ln.ld %>% 
    tidyr::gather(key = "PC", value = "Loading", -chem) %>% 
    tibble::as_tibble() %>% 
    dplyr::right_join(., colgroups, by = "chem")
  plot_loadings_pca$chem <- factor(as.character(plot_loadings_pca$chem), levels = unique(as.character(plot_loadings_pca$chem)))
  
  loadings <- plot_loadings_pca %>%
    dplyr::filter(PC %in% pcs) %>% 
    ggplot(aes(x = chem, y = Loading, color = !!sym(grouping))) + 
    geom_point() +
    geom_segment(aes(yend=0, xend = chem)) +
    facet_wrap(~ PC) + theme_bw() +
    theme(legend.position = "bottom", 
          axis.text.x = element_text(angle = 45, hjust = 1),
          strip.background = element_rect(fill = "white"),
          axis.title.x = element_blank(),
          axis.title.y = element_blank()) +
    geom_hline(yintercept = 0, linewidth = 0.2) + 
    ggtitle("Principle Component Loadings")
  
  list(var = var_tbl, var_df = pca_summary, load = loadings)
}

prelim_pca <- pca(queens_standardized)

Now let’s take a look at our preliminary PCA. We’ll start by examining the singular values and the proportion of variance explained by each component:

prelim_pca$var

prelim_pca$var_df %>% 
  ggplot(aes(x = Principle.component, y = Singular.values)) +
  geom_point() +
  geom_segment(aes(y = Singular.values, yend = 0, xend = Principle.component)) +
  theme_bw() +
  labs(x = "Principal component", y = "Singular value") +
  ggtitle("Singular values of preliminary PCA")

And next let’s examine the top few PCA loadings, corresponding to each of the first six principal components (PCs), or first six patterns governing our data:

prelim_pca$load

3.3 PCP parameter gridsearch

Now we’re about ready to run our first gridsearch.

D <- queens_standardized %>% 
  select(-Date) %>% 
  as.matrix() %>% 
  apply(., 2, mean_impute)

We can run a grid search using the default_eta and searching through all ranks 1 through 15 with the following code chunk:

pcp_defaults <- queens_standardized %>% select(-Date) %>% get_pcp_defaults()
default_eta <- pcp_defaults$eta

grid <- grid_search_cv(
  mat = D, 
  pcp_func = RRMC, 
  grid = data.frame(r = c(15)), 
  eta = default_eta
)

This takes quite a while, so we have saved the results to a .rds file using grid_search_cv’s save_as argument, which we load below:

grid <- readRDS(here::here("my-doc", "huge_eta_search.rds"))
grid$all_stats %>% 
  group_by(r, eta) %>% 
  summarize(
    avg_rel_err = mean(rel_err), 
    avg_sparsity = mean(S_sparsity), 
    avg_rank = mean(L_rank)
  ) %>% 
  arrange(avg_rel_err)

## # A tibble: 66 × 5
## # Groups:   r [6]
##        r   eta avg_rel_err avg_sparsity avg_rank
##    <int> <dbl>       <dbl>        <dbl>    <dbl>
##  1     3  0.2        0.846        100.         3
##  2     3  0.15       0.846         99.9        3
##  3     3  0.1        0.852         99.8        3
##  4     3  0.3        0.853        100.         3
##  5     3  0.25       0.856        100.         3
##  6     2  0.3        0.881        100.         2
##  7     3  0.4        0.882        100.         3
##  8     2  0.25       0.882        100.         2
##  9     3  0.35       0.882        100.         3
## 10     2  0.2        0.884        100.         2
## # ℹ 56 more rows

plot_mat <- function(D) {
  heatmaply::heatmaply(D, Rowv=F, Colv=F, showticklabels = FALSE, showtickmarks = FALSE)
}

rrmc_out <- RRMC(D, r = 3, eta = 0.1)

plot_mat(rrmc_out$L)

plot_mat(rrmc_out$S)

colnames(rrmc_out$L) <- colnames(D)
rrmc_pca <- pca(rrmc_out$L, pcs = paste("PC", 1:3, sep = ""))
rrmc_pca$var

rrmc_pca$load

print_patterns <- function(pats, colgroups = NULL, n = 1:6, pat_type = "pat", title = "") {
  
  if (!is.null(colgroups)) {
    colgroups <- colgroups %>% dplyr::rename(chem = !!names(colgroups)[1])
  } else {
    colgroups <- data.frame(chem = rownames(pats), group = "1")
  }
  
  grouping <- names(colgroups)[2]
  
  colnames(pats) <- paste0(pat_type, stringr::str_pad(1:ncol(pats), width = 2, pad = "0", side = "left"))
  
  pats.df <- pats %>% 
    tibble::as_tibble() %>% 
    dplyr::mutate(chem = colgroups[[1]]) %>%
    tidyr::pivot_longer(-chem, names_to = "pattern", values_to = "loading") %>%
    dplyr::right_join(., colgroups, by = "chem")

  pats.df$chem <- factor(as.character(pats.df$chem), levels = unique(as.character(pats.df$chem)))

  loadings <- pats.df %>%
    dplyr::filter(pattern %in% paste0(pat_type, stringr::str_pad(n, width = 2, pad = "0", side = "left"))) %>%
    ggplot(aes(x = chem, y = loading, color = !!sym(grouping))) +
    geom_point() +
    geom_segment(aes(yend=0, xend = chem)) +
    facet_wrap(~ pattern) + theme_bw() +
    theme(legend.position = "bottom",
          axis.text.x = element_text(angle = 45, hjust = 1),
          strip.background = element_rect(fill = "white"),
          axis.title.x = element_blank(),
          axis.title.y = element_blank()) +
    geom_hline(yintercept = 0, size = 0.2) + ggtitle(title)
  
  loadings
}

# for FA
library(psych)
library(GPArotation)
library(ggrepel)

rrmc_fa <- fa(rrmc_out$L, nfactors = 3, n.obs = nrow(rrmc_out$L), rotate = "varimax", scores = "regression")

## Warning in cor.smooth(r): Matrix was not positive definite, smoothing was done

print(rrmc_fa, digits = 2)

## Factor Analysis using method =  minres
## Call: fa(r = rrmc_out$L, nfactors = 3, n.obs = nrow(rrmc_out$L), rotate = "varimax", 
##     scores = "regression")
## Standardized loadings (pattern matrix) based upon correlation matrix
##       MR1   MR2   MR3 h2      u2 com
## Al  -0.05  0.98  0.20  1 0.00076 1.1
## NH4  1.00  0.09 -0.01  1 0.00036 1.0
## As   1.00  0.04  0.02  1 0.00038 1.0
## Ba  -0.18  0.97 -0.14  1 0.00086 1.1
## Br   0.95  0.16  0.27  1 0.00043 1.2
## Cd   0.83 -0.13 -0.55  1 0.00067 1.8
## Ca   0.77  0.60  0.23  1 0.00038 2.1
## Cl   0.10 -0.07  0.99  1 0.00137 1.0
## Cr   0.38  0.92  0.12  1 0.00056 1.4
## Cu   0.69  0.71 -0.16  1 0.00042 2.1
## EC   0.83  0.53 -0.17  1 0.00035 1.8
## Fe   0.56  0.82 -0.06  1 0.00047 1.8
## Pb   0.97  0.04 -0.26  1 0.00042 1.1
## Mg  -0.12  0.17  0.98  1 0.00128 1.1
## Mn   0.77  0.64  0.04  1 0.00035 1.9
## Ni   1.00  0.03  0.06  1 0.00039 1.0
## OC   0.77  0.63 -0.14  1 0.00037 2.0
## K    0.90  0.34  0.25  1 0.00038 1.5
## Se   0.97  0.00 -0.22  1 0.00042 1.1
## Si   0.06  0.98  0.16  1 0.00072 1.1
## Na   0.06  0.02  1.00  1 0.00135 1.0
## S    0.95  0.30  0.07  1 0.00032 1.2
## Ti  -0.06  1.00 -0.02  1 0.00079 1.0
## NO3  1.00 -0.03  0.05  1 0.00042 1.0
## V    0.99  0.10  0.10  1 0.00038 1.0
## Zn   1.00  0.02 -0.03  1 0.00039 1.0
## 
##                         MR1  MR2  MR3
## SS loadings           14.62 7.63 3.74
## Proportion Var         0.56 0.29 0.14
## Cumulative Var         0.56 0.86 1.00
## Proportion Explained   0.56 0.29 0.14
## Cumulative Proportion  0.56 0.86 1.00
## 
## Mean item complexity =  1.3
## Test of the hypothesis that 3 factors are sufficient.
## 
## df null model =  325  with the objective function =  523.71 with Chi Square =  798399
## df of  the model are 250  and the objective function was  331.99 
## 
## The root mean square of the residuals (RMSR) is  0 
## The df corrected root mean square of the residuals is  0 
## 
## The harmonic n.obs is  1535 with the empirical chi square  0.1  with prob <  1 
## The total n.obs was  1535  with Likelihood Chi Square =  505460.2  with prob <  0 
## 
## Tucker Lewis Index of factoring reliability =  0.176
## RMSEA index =  1.147  and the 90 % confidence intervals are  1.145 1.15
## BIC =  503626.1
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy             
##                                                   MR1 MR2 MR3
## Correlation of (regression) scores with factors     1   1   1
## Multiple R square of scores with factors            1   1   1
## Minimum correlation of possible factor scores       1   1   1

loadings <- as.data.frame(cbind(rownames(rrmc_fa$loadings[]), rrmc_fa$loadings[])) %>% rename(Variable = V1) 

loadings <- loadings %>% mutate_at(colnames(loadings)[str_starts(colnames(loadings), "MR")], as.numeric)

loadings$Max <- colnames(loadings[, -1])[max.col(loadings[, -1], ties.method = "first")] # should be 2:5

loadings %>% kableExtra::kbl(caption = "Loadings") %>% kableExtra::kable_classic(full_width = F, html_font = "Cambria", position = "center") %>% 
  kableExtra::kable_styling(bootstrap_options = c("hover", "condensed"), fixed_thead = T) %>% kableExtra::scroll_box(width = "100%", height = "400px")

scores <- as.data.frame(cbind(rownames(rrmc_fa$scores[]), rrmc_fa$scores[])) %>% mutate_all(as.numeric)

scores$Max <- colnames(scores)[max.col(scores, ties.method = "first")]

scores %>% 
  kableExtra::kbl(caption = "Scores") %>% 
  kableExtra::kable_classic(full_width = F, html_font = "Cambria", position = "center") %>% 
  kableExtra::kable_styling(bootstrap_options = c("hover", "condensed"), fixed_thead = T) %>% 
  kableExtra::scroll_box(width = "100%", height = "400px")